#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Feb 17 11:44:57 2024

@author: jcfq2
"""

import os

jwst_dir='/Users/jcfq2/data/observations/jwst'

os.chdir(jwst_dir)

import matplotlib.pyplot as plt
import matplotlib.colors as colours
from astropy.visualization import make_lupton_rgb
from astropy.table import Table
import ch4_fiddlesticks as ch4
import pandas as pd
import h3ppy
import glob
import spiceypy as spice
from astropy.io import fits
import numpy as np
from importlib import reload
import JWSTSolarSystemPointing as jssp
reload(jssp)
#import jwst_uranus as jwstu

# In[0]

kernel_dir = '/Users/jcfq2/data/observations/jwst/kernels/'
jssp.load_kernels(kdir=kernel_dir)

# Set up a h3ppy object too, always useful
h3p_model = h3ppy.h3p()

# In[1]

fn = 'jw05308001001_03119_g395h-f290lp_s3d.fits'
file_dir = '/Users/jcfq2/data/observations/jwst/5308_dither_combined/'

file = file_dir + fn
# Create a JWSTSolarSystemPointing object that helps with lots of things JWST
geo = jssp.JWSTSolarSystemPointing(file)

# Caluclate the geometry for the full IFU cube. This is a three dimensional arrray.
cube = geo.full_fov()

# Get the wavelength scale for this observation - should be the same for all NIRSpec G395H/F290L observations.
wave = geo.get_wavelength()

# These are the available geometric parameters - if you need something that's not here, let me know!
# Only using Pandas to make this table pretty, generally not a fan
pd.DataFrame(geo.keys, columns=["Parameter"])




# %% test the frame sizes in axes seetting


# plt.figure()

# a3x1_base = plt.subplot(ynum, 5, ypos+(6-i))
# # # gl = a3x1_base.gridlines(linewidth=1, color='grey', alpha=0.0, linestyle='dotted', draw_labels=False)

# # a3x1_chartBox = a3x1_base.get_position()

# # a3x1_img = fig.add_axes([a3x1_chartBox.x0, a3x1_chartBox.y0, 
# #              a3x1_chartBox.width, 
# #              a3x1_chartBox.height], transform=a3x1_base.transAxes, frameon=False) 
# # a3x1_img.set_axis_off()
# # a3x1_base.set_axis_off()

# # a3x1_img.plot(np.arange(100))

# # a3x1 = fig.add_axes([a3x1_chartBox.x0, a3x1_chartBox.y0, 
# #                   a3x1_chartBox.width, 
# #                   a3x1_chartBox.height], transform=a3x1_base.transAxes, frameon=False) 

# a3x1_img = fig.add_axes([0,0,0.5,0.5])
# a3x1_img.plot(np.arange(100),linewidth=2,alpha=0.5)

# # ax1 = fig.add_axes([x0,0, xwidth,  ywidth])   #, frameon=False

# # ax1.plot(np.arange(100))

# plt.show()
# %%



files = sorted(glob.glob("/Users/jcfq2/data/observations/jwst/5308_dither_combined/*.fits"))
geo = jssp.JWSTSolarSystemPointing(files[0])
start1 = geo.obs_start[11:16]
geo = jssp.JWSTSolarSystemPointing(files[19])
start2 = geo.obs_start[11:16]
geo = jssp.JWSTSolarSystemPointing(files[21])
start3 = geo.obs_start[11:16]
geo = jssp.JWSTSolarSystemPointing(files[23])
start4 = geo.obs_start[11:16]
geo = jssp.JWSTSolarSystemPointing(files[25])
start5 = geo.obs_start[11:16]


scale=2
scaleF=1


enc_int=np.load('enceladus_spot_presub.npy')
enc_diff=np.load('enceladus_spot.npy')

nam0=np.load('saturn_spectra_map_1.npy')
nac0=np.load('saturn_spectra_count_1.npy')


nam0[-90*scaleF:-45*scaleF,:,:]=nam0[-90*scaleF:-45*scaleF,:,:]+nam0[0:45*scaleF,:,:]
nac0[-90*scaleF:-45*scaleF,:,:]=nac0[-90*scaleF:-45*scaleF,:,:]+nac0[0:45*scaleF,:,:]
nam0[45*scaleF:90*scaleF,:,:]=nam0[45*scaleF:90*scaleF,:,:]+nam0[-45*scaleF:,:,:]
nac0[45*scaleF:90*scaleF,:,:]=nac0[45*scaleF:90*scaleF,:,:]+nac0[-45*scaleF:,:,:]
ab0=nam0[:,:,2239]+nam0[:,:,1568]+nam0[:,:,1642]+nam0[:,:,1010]+nam0[:,:,1011]+nam0[:,:,1018]+nam0[:,:,1019]
cd0=(nam0[:,:,2241]*0.8*0.5+nam0[:,:,2236]*1.19*0.5)+((nam0[:,:,1571]+nam0[:,:,1565])*0.5)+(0.4*nam0[:,:,1640]/5.21)+(nam0[:,:,1004]+nam0[:,:,1005]*2)
img_h3p0=((ab0-cd0)/nac0[:,:,600])
map_h3p0=(img_h3p0[45*scaleF:405*scaleF,:].T)

maxint=np.nanmax(map_h3p0)

map_h3p0=np.fliplr(img_h3p0[45*scaleF:405*scaleF,:].T)/maxint

map_h3p0=np.roll(map_h3p0,-90*scaleF)

import cartopy.crs as ccrs

crs = ccrs.NorthPolarStereo(globe=ccrs.Globe(flattening=(0.0)))


namA=np.load('saturn_spectra_map_'+str(scale)+'_1.npy')
nacA=np.load('saturn_spectra_count_'+str(scale)+'_1.npy')

namA[-90*scale:-45*scale,:,:]=namA[-90*scale:-45*scale,:,:]+namA[0:45*scale,:,:]
nacA[-90*scale:-45*scale,:,:]=nacA[-90*scale:-45*scale,:,:]+nacA[0:45*scale,:,:]
namA[45*scale:90*scale,:,:]=namA[45*scale:90*scale,:,:]+namA[-45*scale:,:,:]
nacA[45*scale:90*scale,:,:]=nacA[45*scale:90*scale,:,:]+nacA[-45*scale:,:,:]
abA=namA[:,:,2239]+namA[:,:,1568]+namA[:,:,1642]+namA[:,:,1010]+namA[:,:,1011]+namA[:,:,1018]+namA[:,:,1019]
cdA=(namA[:,:,2241]*0.8*0.5+namA[:,:,2236]*1.19*0.5)+((namA[:,:,1571]+namA[:,:,1565])*0.5)+(0.4*namA[:,:,1640]/5.21)+(namA[:,:,1004]+namA[:,:,1005]*2)
img_h3pA=((abA-cdA)/nacA[:,:,600])
map_h3pA=(img_h3pA[45*scale:405*scale,:].T)/maxint

namZ=np.load('saturn_spectra_map_'+str(scale)+'_25.npy')
nacZ=np.load('saturn_spectra_count_'+str(scale)+'_25.npy')

namZ[-90*scale:-45*scale,:,:]=namZ[-90*scale:-45*scale,:,:]+namZ[0:45*scale,:,:]
nacZ[-90*scale:-45*scale,:,:]=nacZ[-90*scale:-45*scale,:,:]+nacZ[0:45*scale,:,:]
namZ[45*scale:90*scale,:,:]=namZ[45*scale:90*scale,:,:]+namZ[-45*scale:,:,:]
nacZ[45*scale:90*scale,:,:]=nacZ[45*scale:90*scale,:,:]+nacZ[-45*scale:,:,:]
abZ=namZ[:,:,2239]+namZ[:,:,1568]+namZ[:,:,1642]+namZ[:,:,1010]+namZ[:,:,1011]+namZ[:,:,1018]+namZ[:,:,1019]
cdZ=(namZ[:,:,2241]*0.8*0.5+namZ[:,:,2236]*1.19*0.5)+((namZ[:,:,1571]+namZ[:,:,1565])*0.5)+(0.4*namZ[:,:,1640]/5.21)+(namZ[:,:,1004]+namZ[:,:,1005]*2)
img_h3pZ=((abZ-cdZ)/nacZ[:,:,600])
map_h3pZ=(img_h3pZ[45*scale:405*scale,:].T)/maxint


fig2, axs = plt.subplots(5, 3, figsize=(11, 7.5), constrained_layout=True,dpi=300)

# fig2 = plt.figure(figsize=(12,6.5),dpi=300)

fig2.subplots_adjust(hspace=0,wspace=0.0)
bbox_props1= dict(boxstyle="circle,pad=0.15", fc="whitesmoke", ec="silver", lw=2)



longmax=360*scale-0*scale
longmin=360*scale-60*scale



ax01 = plt.subplot(5,3,1,projection=ccrs.NorthPolarStereo(central_longitude=180))

ax01.set_extent([100, 270, 40, 90], ccrs.PlateCarree())
ish3p2=ax01.imshow(map_h3p0, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),vmax=0.16,vmin=0.06)
# ax.gridlines()
ax01.text(360-240,60,r'$\alpha$', transform=ccrs.PlateCarree(), ha="center", va="center", weight='bold',c='w',size='small')
ax01.text(360-190,60,r'$\beta$', transform=ccrs.PlateCarree(), ha="center", va="center", weight='bold',c='w',size='small')
ax01.text(360-155,58,r'$\gamma$', transform=ccrs.PlateCarree(), ha="center", va="center", weight='bold',c='w',size='small')
ax01.text(360-132,58,r'$\epsilon$', transform=ccrs.PlateCarree(), ha="center", va="center", weight='bold',c='w',size='small')
# cbar=fig2.colorbar(ish3p2,location='bottom',aspect=8,pad=0.01,shrink=0.85)
# cbar.ax.set_xlabel('H$_3^+$ [mW/m$^2$/sr]')
axs[0,0].set_axis_off()


ax02 = plt.subplot(5,3,2)
ax02.imshow(enc_int/maxint,vmax=0.16,vmin=0.06,extent=[-10,10,60,70],origin='lower')
# ax01.text(longmin/scale,50.5,str(start1)+'UT',va='bottom',ha='right',c='k')

import matplotlib.patches as patches
rect = patches.Rectangle((-2, 64.8-2), 4, 4, linewidth=1, edgecolor='white',facecolor='none')
ax02.add_patch(rect)
ax02.text(2.2,64.8,'Enceladus',va='center',ha='left',c='white')


# ax12 = plt.subplot(5,3,5)
# ax12.imshow(np.fliplr(map_h3pZ[140*scale:155*scale,longmin:longmax]),vmax=0.16,extent=[longmax/scale,longmin/scale,50,65],origin='lower')
# ax12.sharey(ax11)
# ax12.tick_params(left = False,labelleft=False) 
# ax12.text(longmin/scale,50.5,str(start5)+'UT',va='bottom',ha='right')

ax03 = plt.subplot(5,3,3)
ax03.imshow(enc_diff/maxint,vmin=-0.05,vmax=0.05,cmap='seismic',extent=[-10,10,60,70],origin='lower')

ax03.sharey(ax02)
ax03.tick_params(left = False,labelleft=False) 
# ax03.text(longmin/scale,50.5,str(int(start5[0:2])*60-int(start1[0:2])*60+int(start5[3:5])-int(start1[3:5]))+' mins',va='bottom',ha='right',c='k')
rect2 = patches.Rectangle((-2, 64.8-2), 4, 4, linewidth=1, edgecolor='black',facecolor='none')
ax03.add_patch(rect2)
ax03.text(2.2,64.8,'Enceladus',va='center',ha='left',c='black')



diffAZ = map_h3pA[140*scale:155*scale,longmin:longmax]-map_h3pZ[140*scale:155*scale,longmin:longmax]
ax11 = plt.subplot(5,3,4)
ax11.imshow(np.fliplr(map_h3pA[140*scale:155*scale,longmin:longmax]),vmax=0.16,extent=[longmax/scale,longmin/scale,50,65],origin='lower')
ax11.text(longmin/scale,50.5,str(start1)+'UT',va='bottom',ha='right',c='k')

ax11.text(240+90,60,r'$\alpha$', ha="center", va="center", weight='bold',c='w',size='large')
# ax11.set_xticks()

ax12 = plt.subplot(5,3,5)
ax12.imshow(np.fliplr(map_h3pZ[140*scale:155*scale,longmin:longmax]),vmax=0.16,extent=[longmax/scale,longmin/scale,50,65],origin='lower')
ax12.sharey(ax11)
ax12.tick_params(left = False,labelleft=False) 
ax12.text(longmin/scale,50.5,str(start5)+'UT',va='bottom',ha='right',c='w')
ax12.text(240+90,60,r'$\alpha$', ha="center", va="center", weight='bold',c='w',size='large')

ax13 = plt.subplot(5,3,6)
ax13.imshow(np.fliplr(diffAZ),vmin=-0.05,vmax=0.05,cmap='seismic',extent=[longmax/scale,longmin/scale,50,65],origin='lower')

ax13.sharey(ax11)
ax13.tick_params(left = False,labelleft=False) 
ax13.text(longmin/scale,50.5,str(int(start5[0:2])*60-int(start1[0:2])*60+int(start5[3:5])-int(start1[3:5]))+' mins',va='bottom',ha='right',c='k')
ax13.text(240+90,60,r'$\alpha$', ha="center", va="center", weight='bold',c='k',size='large')



namA=np.load('saturn_spectra_map_'+str(scale)+'_23.npy')
nacA=np.load('saturn_spectra_count_'+str(scale)+'_23.npy')

namA[-90*scale:-45*scale,:,:]=namA[-90*scale:-45*scale,:,:]+namA[0:45*scale,:,:]
nacA[-90*scale:-45*scale,:,:]=nacA[-90*scale:-45*scale,:,:]+nacA[0:45*scale,:,:]
namA[45*scale:90*scale,:,:]=namA[45*scale:90*scale,:,:]+namA[-45*scale:,:,:]
nacA[45*scale:90*scale,:,:]=nacA[45*scale:90*scale,:,:]+nacA[-45*scale:,:,:]
abA=namA[:,:,2239]+namA[:,:,1568]+namA[:,:,1642]+namA[:,:,1010]+namA[:,:,1011]+namA[:,:,1018]+namA[:,:,1019]
cdA=(namA[:,:,2241]*0.8*0.5+namA[:,:,2236]*1.19*0.5)+((namA[:,:,1571]+namA[:,:,1565])*0.5)+(0.4*namA[:,:,1640]/5.21)+(namA[:,:,1004]+namA[:,:,1005]*2)
img_h3pA=((abA-cdA)/nacA[:,:,600])
map_h3pA=(img_h3pA[45*scale:405*scale,:].T)/maxint

namZ=np.load('saturn_spectra_map_'+str(scale)+'_25.npy')
nacZ=np.load('saturn_spectra_count_'+str(scale)+'_25.npy')

namZ[-90*scale:-45*scale,:,:]=namZ[-90*scale:-45*scale,:,:]+namZ[0:45*scale,:,:]
nacZ[-90*scale:-45*scale,:,:]=nacZ[-90*scale:-45*scale,:,:]+nacZ[0:45*scale,:,:]
namZ[45*scale:90*scale,:,:]=namZ[45*scale:90*scale,:,:]+namZ[-45*scale:,:,:]
nacZ[45*scale:90*scale,:,:]=nacZ[45*scale:90*scale,:,:]+nacZ[-45*scale:,:,:]
abZ=namZ[:,:,2239]+namZ[:,:,1568]+namZ[:,:,1642]+namZ[:,:,1010]+namZ[:,:,1011]+namZ[:,:,1018]+namZ[:,:,1019]
cdZ=(namZ[:,:,2241]*0.8*0.5+namZ[:,:,2236]*1.19*0.5)+((namZ[:,:,1571]+namZ[:,:,1565])*0.5)+(0.4*namZ[:,:,1640]/5.21)+(namZ[:,:,1004]+namZ[:,:,1005]*2)
img_h3pZ=((abZ-cdZ)/nacZ[:,:,600])
map_h3pZ=(img_h3pZ[45*scale:405*scale,:].T)/maxint


longmax=360*scale-50*scale
longmin=360*scale-110*scale


diffAZ = map_h3pA[140*scale:155*scale,longmin:longmax]-map_h3pZ[140*scale:155*scale,longmin:longmax]
ax21 = plt.subplot(5,3,7)
ax21.imshow(np.fliplr(map_h3pA[140*scale:155*scale,longmin:longmax]),vmax=0.16,extent=[longmax/scale,longmin/scale,50,65],origin='lower')
ax21.text(longmin/scale,50.5,str(start4)+'UT',va='bottom',ha='right',c='w')
ax21.text(190+90,60,r'$\beta$', ha="center", va="center", weight='bold',c='w',size='large')


ax22 = plt.subplot(5,3,8)
ax22.imshow(np.fliplr(map_h3pZ[140*scale:155*scale,longmin:longmax]),vmax=0.16,extent=[longmax/scale,longmin/scale,50,65],origin='lower')
ax22.sharey(ax21)
ax22.tick_params(left = False,labelleft=False) 
ax22.text(longmin/scale,50.5,str(start5)+'UT',va='bottom',ha='right',c='k')
ax22.text(190+90,60,r'$\beta$', ha="center", va="center", weight='bold',c='w',size='large')


ax23 = plt.subplot(5,3,9)
ax23.imshow(np.fliplr(diffAZ),vmin=-0.05,vmax=0.05,cmap='seismic',extent=[longmax/scale,longmin/scale,50,65],origin='lower')
ax23.sharey(ax21)
ax23.tick_params(left = False,labelleft=False) 
ax23.text(longmin/scale,50.5,str(int(start5[0:2])*60-int(start4[0:2])*60+int(start5[3:5])-int(start4[3:5]))+' mins',va='bottom',ha='right',c='k')
ax23.text(190+90,60,r'$\beta$', ha="center", va="center", weight='bold',c='k',size='large')




namA=np.load('saturn_spectra_map_'+str(scale)+'_21.npy')
nacA=np.load('saturn_spectra_count_'+str(scale)+'_21.npy')

namA[-90*scale:-45*scale,:,:]=namA[-90*scale:-45*scale,:,:]+namA[0:45*scale,:,:]
nacA[-90*scale:-45*scale,:,:]=nacA[-90*scale:-45*scale,:,:]+nacA[0:45*scale,:,:]
namA[45*scale:90*scale,:,:]=namA[45*scale:90*scale,:,:]+namA[-45*scale:,:,:]
nacA[45*scale:90*scale,:,:]=nacA[45*scale:90*scale,:,:]+nacA[-45*scale:,:,:]
abA=namA[:,:,2239]+namA[:,:,1568]+namA[:,:,1642]+namA[:,:,1010]+namA[:,:,1011]+namA[:,:,1018]+namA[:,:,1019]
cdA=(namA[:,:,2241]*0.8*0.5+namA[:,:,2236]*1.19*0.5)+((namA[:,:,1571]+namA[:,:,1565])*0.5)+(0.4*namA[:,:,1640]/5.21)+(namA[:,:,1004]+namA[:,:,1005]*2)
img_h3pA=((abA-cdA)/nacA[:,:,600])
map_h3pA=(img_h3pA[45*scale:405*scale,:].T)/maxint

namZ=np.load('saturn_spectra_map_'+str(scale)+'_23.npy')
nacZ=np.load('saturn_spectra_count_'+str(scale)+'_23.npy')

namZ[-90*scale:-45*scale,:,:]=namZ[-90*scale:-45*scale,:,:]+namZ[0:45*scale,:,:]
nacZ[-90*scale:-45*scale,:,:]=nacZ[-90*scale:-45*scale,:,:]+nacZ[0:45*scale,:,:]
namZ[45*scale:90*scale,:,:]=namZ[45*scale:90*scale,:,:]+namZ[-45*scale:,:,:]
nacZ[45*scale:90*scale,:,:]=nacZ[45*scale:90*scale,:,:]+nacZ[-45*scale:,:,:]
abZ=namZ[:,:,2239]+namZ[:,:,1568]+namZ[:,:,1642]+namZ[:,:,1010]+namZ[:,:,1011]+namZ[:,:,1018]+namZ[:,:,1019]
cdZ=(namZ[:,:,2241]*0.8*0.5+namZ[:,:,2236]*1.19*0.5)+((namZ[:,:,1571]+namZ[:,:,1565])*0.5)+(0.4*namZ[:,:,1640]/5.21)+(namZ[:,:,1004]+namZ[:,:,1005]*2)
img_h3pZ=((abZ-cdZ)/nacZ[:,:,600])
map_h3pZ=(img_h3pZ[45*scale:405*scale,:].T)/maxint


longmax=360*scale-70*scale
longmin=360*scale-130*scale

diffAZ = map_h3pA[140*scale:155*scale,longmin:longmax]-map_h3pZ[140*scale:155*scale,longmin:longmax]
ax31 = plt.subplot(5,3,10)
ax31.imshow(np.fliplr(map_h3pA[140*scale:155*scale,longmin:longmax]),vmax=0.16,extent=[longmax/scale,longmin/scale,50,65],origin='lower')
ax31.text(longmin/scale,50.5,str(start3)+'UT',va='bottom',ha='right',c='w')
ax31.text(190+90,60,r'$\beta$', ha="center", va="center", weight='bold',c='w',size='large')
ax31.text(155+90,58,r'$\gamma$', ha="center", va="center", weight='bold',c='w',size='large')

ax32 = plt.subplot(5,3,11)
ax32.imshow(np.fliplr(map_h3pZ[140*scale:155*scale,longmin:longmax]),vmax=0.16,extent=[longmax/scale,longmin/scale,50,65],origin='lower')
ax32.sharey(ax21)
ax32.tick_params(left = False,labelleft=False) 
ax32.text(longmin/scale,50.5,str(start4)+'UT',va='bottom',ha='right',c='k')
ax32.text(190+90,60,r'$\beta$', ha="center", va="center", weight='bold',c='w',size='large')
ax32.text(155+90,58,r'$\gamma$', ha="center", va="center", weight='bold',c='w',size='large')

ax33 = plt.subplot(5,3,12)
ax33.imshow(np.fliplr(diffAZ),vmin=-0.05,vmax=0.05,cmap='seismic',extent=[longmax/scale,longmin/scale,50,65],origin='lower')

ax33.sharey(ax21)
ax33.tick_params(left = False,labelleft=False) 
ax33.text(longmin/scale,50.5,str(int(start4[0:2])*60-int(start3[0:2])*60+int(start4[3:5])-int(start3[3:5]))+' mins',va='bottom',ha='right',c='k')
ax33.text(190+90,60,r'$\beta$', ha="center", va="center", weight='bold',c='k',size='large')
ax33.text(155+90,58,r'$\gamma$', ha="center", va="center", weight='bold',c='k',size='large')



namA=np.load('saturn_spectra_map_'+str(scale)+'_19.npy')
nacA=np.load('saturn_spectra_count_'+str(scale)+'_19.npy')

namA[-90*scale:-45*scale,:,:]=namA[-90*scale:-45*scale,:,:]+namA[0:45*scale,:,:]
nacA[-90*scale:-45*scale,:,:]=nacA[-90*scale:-45*scale,:,:]+nacA[0:45*scale,:,:]
namA[45*scale:90*scale,:,:]=namA[45*scale:90*scale,:,:]+namA[-45*scale:,:,:]
nacA[45*scale:90*scale,:,:]=nacA[45*scale:90*scale,:,:]+nacA[-45*scale:,:,:]
abA=namA[:,:,2239]+namA[:,:,1568]+namA[:,:,1642]+namA[:,:,1010]+namA[:,:,1011]+namA[:,:,1018]+namA[:,:,1019]
cdA=(namA[:,:,2241]*0.8*0.5+namA[:,:,2236]*1.19*0.5)+((namA[:,:,1571]+namA[:,:,1565])*0.5)+(0.4*namA[:,:,1640]/5.21)+(namA[:,:,1004]+namA[:,:,1005]*2)
img_h3pA=((abA-cdA)/nacA[:,:,600])
map_h3pA=(img_h3pA[45*scale:405*scale,:].T)/maxint

namZ=np.load('saturn_spectra_map_'+str(scale)+'_21.npy')
nacZ=np.load('saturn_spectra_count_'+str(scale)+'_21.npy')

namZ[-90*scale:-45*scale,:,:]=namZ[-90*scale:-45*scale,:,:]+namZ[0:45*scale,:,:]
nacZ[-90*scale:-45*scale,:,:]=nacZ[-90*scale:-45*scale,:,:]+nacZ[0:45*scale,:,:]
namZ[45*scale:90*scale,:,:]=namZ[45*scale:90*scale,:,:]+namZ[-45*scale:,:,:]
nacZ[45*scale:90*scale,:,:]=nacZ[45*scale:90*scale,:,:]+nacZ[-45*scale:,:,:]
abZ=namZ[:,:,2239]+namZ[:,:,1568]+namZ[:,:,1642]+namZ[:,:,1010]+namZ[:,:,1011]+namZ[:,:,1018]+namZ[:,:,1019]
cdZ=(namZ[:,:,2241]*0.8*0.5+namZ[:,:,2236]*1.19*0.5)+((namZ[:,:,1571]+namZ[:,:,1565])*0.5)+(0.4*namZ[:,:,1640]/5.21)+(namZ[:,:,1004]+namZ[:,:,1005]*2)
img_h3pZ=((abZ-cdZ)/nacZ[:,:,600])
map_h3pZ=(img_h3pZ[45*scale:405*scale,:].T)/maxint


longmax=360*scale-90*scale
longmin=360*scale-150*scale

diffAZ = map_h3pA[140*scale:155*scale,longmin:longmax]-map_h3pZ[140*scale:155*scale,longmin:longmax]
ax41 = plt.subplot(5,3,13)
imh3pA=ax41.imshow(np.fliplr(map_h3pA[140*scale:155*scale,longmin:longmax]),vmax=0.16,extent=[longmax/scale,longmin/scale,50,65],origin='lower')
ax41.text(longmin/scale,50.5,str(start2)+'UT',va='bottom',ha='right',c='w')
ax41.text(155+90,58,r'$\gamma$', ha="center", va="center", weight='bold',c='w',size='large')
ax41.text(132+90,58,r'$\epsilon$', ha="center", va="center", weight='bold',c='w',size='large')

ax42 = plt.subplot(5,3,14)
imh3pZ=ax42.imshow(np.fliplr(map_h3pZ[140*scale:155*scale,longmin:longmax]),vmax=0.16,extent=[longmax/scale,longmin/scale,50,65],origin='lower')

ax42.sharey(ax21)
ax42.tick_params(left = False,labelleft=False) 
ax42.text(longmin/scale,50.5,str(start3)+'UT',va='bottom',ha='right',c='w')
ax42.text(155+90,58,r'$\gamma$', ha="center", va="center", weight='bold',c='w',size='large')
ax42.text(132+90,58,r'$\epsilon$', ha="center", va="center", weight='bold',c='w',size='large')

ax43 = plt.subplot(5,3,15)
imAZ=ax43.imshow(np.fliplr(diffAZ),vmin=-0.05,vmax=0.05,cmap='seismic',extent=[longmax/scale,longmin/scale,50,65],origin='lower')

ax43.sharey(ax21)
ax43.tick_params(left = False,labelleft=False) 
ax43.text(longmin/scale,50.5,str(int(start3[0:2])*60-int(start2[0:2])*60+int(start3[3:5])-int(start2[3:5]))+' mins',va='bottom',ha='right',c='k')
ax43.text(155+90,58,r'$\gamma$', ha="center", va="center", weight='bold',c='k',size='large')
ax43.text(132+90,58,r'$\epsilon$', ha="center", va="center", weight='bold',c='k',size='large')


cbar=fig2.colorbar(imh3pA, ax=axs[:, 0:2], shrink=0.8,location='bottom')
cbar.ax.set_xlabel('H$_3^+$')

cbar=fig2.colorbar(imAZ, ax=axs[:, -1], shrink=0.8,location='bottom',aspect=7.6)
cbar.ax.set_xlabel('diff H$_3^+$')


fig2.savefig('asd_fig3.pdf', dpi=300, bbox_inches='tight', facecolor='white', pad_inches=0) 



